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^  ABSTRACT 

The  ability  to  model  random  time  aarlaa  plays 
a  proalncnt  rola  In  a  variety  of  applications  as 
exsapllfled  by  aalamlc  data  analysis,  dopplar 
radar  processing,  speech  processing,  adaptive 
filtering,  and,  array  procasslng.  Undoubtedly, 
two  of  the  more  popular  procedures  for  effecting 
such  time  series  models  are  the  classical  Fourier 
(MA)  approach  and  the  maxlmun  entropy  (AR)  method. 
In  this  paper,  a  theoretical  comparison  of  these 
contemporary  procedures  with  a  more  general  ARKA 
method  will  be  made.  It  will  be  demonstrated  that 
the  spectral  estimation  performance  of  the  ARMA 
method  typically  exceeds  that  of  Its  more  special¬ 
ized  MA  and  AR  counterparts.  With  this  supremacy 
thus  established,  a  recently  developed  method  for 
effectively  generating  ARMA  model  estimates  from 
time  series  observations  will  be  then  presented. 

\ 

I .  INTRODUCTION 

A  problem  which  arises  in  a  variety  of  appli¬ 
cations  Is  that  of  estimating  the  statistical 
characteristics  of  a  random  wlde-sense  stationary 
time  series  tx(n)).  This  estimation  la  typically 
based  upon  a  finite  set  of  time  series  observations, 
in  many  signal  processing  applications,  only  the 
second  order  statistics  as  represented  by  the  time 
series'  autocorrelation  sequence 

r^(n)  •  Eix(n+i!i)x*(m) )  (1) 

is  required.  In  this  expression,  the  symbols  E 
and  *  denote  the  operations  of  expectation  and  com¬ 
plex  conjugation,  respectively.  Upon  taking  the 
Fourier  transform  of  this  deterministic  auto¬ 
correlation  sequence,  we  obtain  the  associated 
power  spectral  density  function 

S  (e-^“)  -  I  r  (n)e  (2) 


It  often  happens  that  the  essential  attributes  of 
a  time  series  are  more  discernible  from  its  fre¬ 
quency  domain  spectral  density  function  than  from 
its  equivalent  time  domain  autocorrelation  sequence. 
It  is  with  this  in  mind  chat  interest  in  spectral 
estimation  techniques  has  evolved. 

This  work  was  supported  in  part  by  the  Office  of 
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In  the  classical  estimation  problem.  It  Is 
desired  to  achieve  an  estimation  of  the  spectral 
density  function  (2)  from  observations  of  a  time 
series.  Without  loss  of  generality,  these  obser¬ 
vations  will  be  here  taken  to  be  the  N  contiguous 
elements 

x(l),  x(2) . x(N)  (3) 

A  variety  of  procedures  have  been  proposed  for 
using  these  observations  to  achieve  a  spectral 
density  estimate.  Without  doubt,  the  overwhelislng 
number  of  procedures  ultimately  reault  In  a 
rational  spectral  density  model  which  fits  the  form 


I  ,  •  1  (a1 

|l+a^e  + 


A  (e>) 


The  S)^  and  bj^  coefficients  of  this  model  are 
referred  to  as  its  autoregressive  and  moving 
average  coefficients,  respectively.  This  model  Is 
commonly  referred  to  as  an  autoregresslvc-movlng 
average  (ARMA)  spectral  model  of  order  (p,q)  where 
q  and  p  denote  Che  orders  of  the  numerator  and 
denominator  polynomials,  respectively.  It  is 
readily  shown  Chet  any  continuous  (In  u)  spectral 
density  may  be  approximated  arbitrarily  closely  by 
the  above  rational  model  If  the  order  (p,q)  Is 
salacced  suitably  large.  Thus,  the  robustness  of 
this  rational  model  is  apparent. 

In  studies  directed  cowards  spectral  analysis, 

Che  preponderance  of  effort  has  been  directed  to¬ 
wards  two  special  cases  of  the  general  ARMA  model 
(A).  They  are  the  moving  average  (MA)  model  for  _ 
which  Ap(eJ")  3  1,  and,  Che  autoregressive  (AR) 
model  for  which  Bq(e)'^).  ■  bg.  The  spectral  density  — 
arising  from  a  HA  model  Is  seen  to  contain  no  poles, 
and,  as  such  it  Is  known  as  an  all-zero  model. 
Similarly,  the  AR  model  Is  referred  to  as  an  all¬ 
pole  model,  and,  the  general  ARMA  model  Is  seen  to 
be  a  pole-zero  model.  Undoubtedly,  the  primary  n. 

reasons  for  Interest  in  the  special  case  MA  and  AR  _ 

models  are  chat  they: 


(1)  are  amenable  to  a  tractable  analysis. 
(11)  typically  provide  adequate  spectral 
estimation  performance. 
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(Ill)  are  syntheslzable  by  efficient  algorlchaa. 

Despite  these  factors,  It  Is  widely  rscognlzsd  Ct>sC 
the  more  robust  ARMA  sodel  typically  provldas  a 
ouch  superior  spectral  estloatlon  perforaance  and 
uses  fewer  oodel  paraoeters  In  the  estimate.  The 
main  impediment  to  Its  use  on  a  wider  scale  has  been 
Che  lack  of  a  specific  procedure  for  obtaining  the 
ARMA  model  parameters  In  a  computationally  effici¬ 
ent  manner.  Recently,  an  ASHA  spectral  modeling 
method  which  possesses  this  computational  capabili¬ 
ty  has  been  developed  (  1  ]  -  (  4 ] .  The  main 
features  of  this  new  procedure  are  outlined  in 
Sections  IV  and  VI. 

In  the  next  three  sections,  theoretical  pro¬ 
cedures  for  generating  HA,  AS,  and  ARMA  models  from 
a  finite  set  of  actual  autocorrelation  values  are 
presented.  This  Is  then  followed  by  an  example 
section  which  treats  the  classical  problem  of 
spectral  estimating  a  time  series  composed  of  two 
sinusoids  In  additive  white  noise.  It  is  there 
shown  that  the  more  general  ARHA  modeling  procedure 
easily  outperforms  Its  special  case  MA  and  AR 
modeling  procedures.  With  the  ARHA  spectral  imidallng 
methods  superior  performance  thereby  demonstrated 
for  this  Idealistic  situation  (l.e.,  actual  auto¬ 
correlation  values  are  given) ,  we  next  direct  our 
attention  to  the  more  practical  problem  of  generat¬ 
ing  ARHA  spectral  estimates  from  a  finite  set  of 
observations  (3).  In  particular,  the  recently 
developed  high  performance  ARHA  modeling  method  Is 
outlined  (2]&[4].  This  novel  procedure  has 
been  found  to  outperform  such  alternate  aRHA 
modeling  procedures  as  the  Box-Jenklns  method  (101 
and  whitening  filter  approaches  (11)4[121. 

II.  MA  SPECTRAL  MODELING 

Thera  exist  a  variety  of  procedures  for  obtain¬ 
ing  a  MA  model  of  a  wide  sense  stationary  time 
series.  These  Include  the  perlodogram  and  differ¬ 
ent  versions  of  the  autocorrelation  method.  In 
each  case,  the  spectral  estimate  will  be  of  the 
form 

S^(eJ“)  -  |bo+b^e‘J“+  .  .  .+b^e'J‘>“|^  (5) 

This  MA  spectral  model  requires  a  rather  large 
value  for  the  order  parameter  q  to  enable  It  to 
achieve  a  desirable  frequency  resolution  perform¬ 
ance.  Unfortunately,  this  requirement  can  lead  to 
a  rather  poor  spectral  estimation  behavior  In  the 
case  of  moderate  length  time  aeries  observations. 

This  behavior  will  be  Illustrated  In  the  numerical 
example  section. 


In  which  w(n)  Is  a  symawtrlc  window  function  which 
Is  selected  to  effect  aosw  desired  behavior.  In 
the  pure  truncated  case,  the  rectangular  window  la 
used  In  which  case  w(n)  ■  1  for  |q|<n.  Since  the 
autocorrelation  sequence  is  a  cosq>lex  conjugate 
synsetrlc  function  of  a  (l.e.,  rx(~>i)  "  r;|*(n)), 
one  can  readily  show  that  expression  (6)  can  be 
equivalently  represented  In  the  form  specified  by 
expression  (S)  whereby  the  prevailing  coefficients 
are  related  by 

q 

w(n)r  (n)  -  I  0<nsq  (7) 

k-n 

Given  the  q-fl  values  of  Che  product  w(n)rj,(n),  one 
may  readily  solve  this  nonlinear  system  of  q  4-  1 
equations  to  obtain  Che  MA  spectral  models  ba  co¬ 
efficients  as  used  In  expression  (S).  Expressions 
(6)  and  (7)  then  consclcuca  a  syscamatlc  procedure 
for  generating  a  MA  modal  of  time  series  baaed  on 
given  autocorrelation  values. 

In  most  applications,  one  has  avallabla  only 
a  sat  of  tlma  series  observations  (3)  (and  not 
autocorrelation  values)  iq>on  which  to  generate  a 
MA  spectral  model.  If  expression  (6)  Is  to  be  used 
for  this  objective.  It  will  than  be  necessary  to 
obtain  estimates  of  the  autocorrelation  elements 
from  the  given  time  series  observations.  The  un¬ 
biased  estimator  as  specified  by 

,  M-n 

rx(n)  •  I  x(n+k)xTk)  0<n<q  (8) 


is  often  used  for  this  purpose  in  which  It  Is 
assvned  chat  the  MA  model  order  Is  such  that  q<N. 
Alternatively,  the  perlodogram  method  has  served 
as  a  useful  procedure  for  effecting  a  MA  spectral 
estimate  from  a  set  of  contiguous  time  series  obsar- 
vgclons  [  6]  ,  The  perlodogram  possesses  the 
additional  advantage  of  being  efficiently  Imple¬ 
mented  by  the  fast  Fourier  transform  algorithm. 

The  Inherent  order  of  a  MA  perlodogram  model  Is  N-1. 

III.  AR  SPECTRAL  MODELING 


In  this  section,  the  method  of  linear  pre¬ 
diction  Is  used  for  generating  an  AR  spectral  model 
associated  with  a  given  time  series.  In  particular, 
the  coefficients  of  the  p'^*'  order  AR  spectral  model 
as  specified  l^y 


or 

O 

l*a^e--’“  +  .  . 

.  +a  e-->P“ 
P 

<9) 


In  this  section,  we  will  be  concerned  with 
evolving  a  MA  spectral  esclswce  procedure  for  the 
special  case  In  which  one  has  available  Che  q-t-l 
autocorrelation  elasiencs  t^(0) ,  r^^d),.  .  ..tj^lq). 

With  this  Information  provided,  a  standard  esti¬ 
mation  Is  generated  by  the  truncated  series  [SI. 

S  (eJ“)  -  I  w(n)r  (nje'-J""  (6) 

*  n— q  * 


*LL1  be  determined  by  solving  a  system  of  p+1 
linear  aquations  In  the  p41  coefficient  unknowns 
•l»  ■2»  •  •  ••  Apt  bQ.  Those  aquations  are  obtain¬ 
ed  by  considering  the  specific  problem  of  pre¬ 
dicting  the  tlma  aeries  element  x(n)  by  a  linear 
combination  of  the  p  most  recant  time  series 

elements  x(n-l) ,x{n-2) . x(n-p).  It  will 

turn  out  chat  the  resultant  sat  of  optimal  equat¬ 
ions  thus  obtained  will  correspond  exactly  with 
those  equations  which  arise  when  using  the 


2 


"maximum  entropy"  method  of  spectral  estimation. 
This  equivalency  has  been  previously  recognized  [  7 ]. 

As  Indicated,  the  task  at  hand  Is  that  of 
estimating  the  time  series  element  x(n)  by  means  of 
the  linear  combination 

P 

x(n)  •  -  I  aj^x(n-k)  (10) 

k»i 

in  which  the  optimal  prediction  coefficients  at^* 
are  to  be  ultimately  used  In  Che  AR  spectral  model 
(9) .  The  prediction  error  is  formally  given  by 


e(n)  -  x(n)  -  x(n) 

P 

•  x(n)  +  I  a^x(n-k)  (11) 

k"l 

It  is  now  desired  to  select  Che  complex  valued  a)^ 
prediction  coefficients  so  chat  the  expected  value 
of  the  error's  magnitude  squared  is  minimized.  One 
may  straightforwardly  show  chat  this  mean  squared 
error  is  given  by 

,  P 

E(,e(n);'='  -r^(0)  +  J^(a^r^(-k)  ^a^r^(k)l 
P  p 

^  J. 

k»i  m”l 

Upon  using  standard  calculus  methods,  it  is  found 
that  the  minimizing  predictor  coefficients  must 
satisfy  the  following  system  of  p  linear  equations 

P 

1  a,"  r  (m-k)  -  -r  (m)  for  l<m<p  (13) 
k-l  “ 


If  this  optimal  selection  is  inserted  into  ex¬ 
pression  (12),  Che  minimum  mean  squared  error  is 
found  to  be 

,  P 

El  a'ln)!^)  -r,(0)  +  la^'r^(-k)  (14) 

X  K  X 


In  suanary,  Che  optimal  predictor  coefficients 
are  obtained  upon  solving  the  system  of  linear 
equations  (13),  and,  Che  corresponding  minimum 
mean  squared  error  is  computed  by  means  of  relation- 
ship  (14).  From  a  computational  viewpoint,  however, 
a  more  efficient  method  for  iteratively  obtaining 
the  optimal  predictor  performance  is  obtained  by 
incorporating  relationships  (13)  and  (14)  into  the 
single  expression 


^r^(O)  r^(-l) 

r  (1)  r  (0) 

’  X  X 


r^(-p) 

r^(-jH.l) 


^1 


(15) 


r^(p) 


r^(p-l)....  r^(0) 


0 


in  which  Ep  -  E(|e'(n)l^}.  One  may  solve  this  Toepllti 
system  of  equations  using  the  Levinson  algorithm 
in  a  computationally  efficient  order  update  manner 


(e.g.,  see  ref.  [7]).  Namely,  Che  optimal  co¬ 
efficients  of  the  ps-l^f  order  predictor  may  be 
recursively  obtained  from  Che  optimal  coefficients 
of  the  pfh  order  predictor.  As  indicated  previous¬ 
ly,  Che  system  of  equations  (IS)  is  identical  Co 
chat  obtained  when  using  Che  maximum  entropy 
method  of  spectral  estimation. 

If  the  optimal  predictor  is  performing  its 
objective.  It  follows  that  Che  prediction  eleawnt 
x(n)  will  contain  all  which  Is  predictable  In  x(n) . 
As  such,  the  prediction  error  (11)  is  white  like 
in  behavior  and  its  spectral  density  is  Chen  given 
by  Se(eJ“)  •  Ep.  This  behavior  la  of  course  de¬ 
pendent  on  making  the  order  parameter  p  sufficient¬ 
ly  large  so  as  to  achieve  Che  desired  perfect 
prediction.  Assuming  this  prediction  behavior,  it 
then  follows  from  relationship  (11)  chat  Che 
spectral  density  function  of  the  time  series  (x(n)) 
is  given  by  expression  (9)  in  which  the  coefficient 
bo^  ■  Ep  and  the  a^^  coefficients  are  obtained 
upon  solving  relationship  (13)  or  equivalently 
relationship  (15). 

In  this  analysis,  it  has  been  assumed  chat 
one  has  available  the  time  series  autocorrelation 
values  rjj(O),  r^d) ,  .  .  .,  rx(p).  More  realistic¬ 
ally,  such  perfect  autocorrelation  knowledge  is 
almost  never  available.  In  a  typical  application, 
one  has  available  only  a  sampled  set  of  time  series 
observations  as  exemplified  by  expression  (3).  If 
the  AR  spectral  estimation  procedure  as  represented 
by  expression  (15)  is  to  be  incorporated,  one  could 
use  Che  given  time  series  observations  to  obtain 
estimates  of  Che  required  autocorrelation  elements. 
Alternatively,  a  set  of  determin< "tic  prediction 
error  equations  can  be  minimized  so  as  to  obtain 
Che  prediction  coefficients.  In  effect,  this  is 
the  approach  usually  taken  in  evolving  the  Burg 
algorithm  [ ^  1 . 


IV.  AJtMA  SPECTRAL  MODELING 

The  time  series  {x(n))  Is  said  to  be  an  ARMA 
process  of  order  (p,q)  If  it  la  generated  according 
to  Che  linear  causal  relationship 

p  q 

x(n)  +  [  a.x(n-k)  •  J  b.w(n-k)  (16) 

k-1  k-0 

where  the  excitation  sequence  {w(n) }  is  a  zero  mean 
white  noise  time  series  whose  Individual  elements 
have  variance  one.  It  Is  a  simple  matter  to  show 
that  the  spectral  density  function  associated  with 
this  response  time  series  (x(n))  Is  given  precisely 
by  expression  (4).  Thus,  there  is  an  equivalency 
between  a  rational  spectral  density  model  and  the 
response  of  a  linear  system  to  a  white  noise 
excitation. 

a^  Coefficient  Determination 

A  procedure  for  identifying  an  ARMA  model's 
a|j  sutoragraaslva  coefficients  Involves  examining 
Its  second  order  statistical  characterization. 

This  is  achieved  by  first  multiplying  both  sides  of 


3 


exprasston  (16)  by  x*(n-»)  and  then  taking  expect¬ 
ed  values.  The  results  of  this  operation  are  the 
well  knotm  Yule-Walker  equations 

P 

r  (m)  +  I  a.  r  (m-k)  “  0  for  m  >  q  (17) 

*  lB-1 

where  It  is  loportant  to  note  that  the  lag  variable 
■  Is  here  restricted  to  be  larger  than  the  ASHA 
Bodel's  nuaerator  order  q. ^ 

In  effect,  the  above  Yule-Walker  equations 
Indicate  that  an  ASHA  time  series'  autocorrelation 
elMMnts  are  Interrelated  in  a  well  defined  linear 
Banner  for  appropriate  lag  values.  This  obser¬ 
vation  than  provides  the  vehicle  for  determining 
the  aSM(  idel's  associated  aj^  coefficients.  To 
be  more  svaclflc,  let  us  now  express  the  first  t 
Yule-Walker  equations  (17)  In  the  following  isatrlx 
format 


t^Cq) 

rx(q-l).  . 

•  •''x^'’+1-p) 

“l 

rx(q+l) 

'rx(q+l) 

rx(q)  .  . 

^2 

m  • 

rx(q+2) 

“p 

rx(<l+'~^) 

rx(q+t-2) . 

•  •'^x^9+t-p) 

r  (q+t) 

* 

(18) 

or  In  the  mors  compact  representation 

R  a  -  -r  (19) 

where  R  la  a  txp  matrix  and  a^  and  £  are  pxl  and 
t*l  vectors,  respectively.  ~ 

It  Is  readily  shown  that  this  system  of  aquations 
will  have  a  unique  solution  provided  that  t  ^  p 
To  obtain  the  ASHA  model's  a^  coefficients,  one 
Chen  simply  solves  this  system  of  linear  equations. 
From  a  computational  viewpoint,  It  Is  convenient 
CO  set  c  equal  to  Its  bIoIbub  value  of  p.  For 
reasons  which  will  be  spelled  out  In  Section  VI, 
however.  It  will  often  be  advantageous  to  let 
c  cake  on  values  exceeding  p. 

b^  Coefficient  Determination 

To  determine  Che  ARMA  model's  b)(  coeffici¬ 
ents,  It  will  be  expedient  to  introduce  the  auto¬ 
correlation  sequence's  causal  Image 

r*(a)  •  rjj(n)u(n)  -  ^  rjj(0)«(n)  (20) 

where  u(n)  and  £(n)  dsnota  tha  standard  unit- 
step  and  Kroneckar  dales  sequencas,  respectively. 
Tha  aucocorralatlon  sequance  may  be  recovered  from 
Its  causal  Image  through  the  relationship 

rx(n)  - 

^Tha  Yule-Walkar  aquations  associated  with  lag 
valuas  O^m^q  will  Involve  tha  ARHA  model's  b)^ 
coafflclMC^  In  a  nonlinear  manner. 


whose  validity  Is  established  by  making  usa  of  tha 
complex  conjugate  synewtry  property  of  auto¬ 
correlation  sequences  (i.e.,  r(-n)  -  r*(n)).  The 
Fourier  transform  of  expression  (21)  Is  found  to 

be 

Sx(e-*“)  -  S^Ce^")  +  Sx^(e->")  (22) 

in  which  5^(0^"’)  denotes  Che  Fourier  transform  of 
Che  autocorrelations'  causal  Image  (20).  In  what 
Is  to  follow,  a  syscamatlc  procedure  for  Identify¬ 
ing  S±(a)“)  will  be  given  which,  with  Che  uclllzsc- 
xon  of  expression  (22),  results  In  Che  ovarall  times 
series'  spectral  density. 


This  Identification  Is  perhaps  best  achieved 
by  Introducing  the  following  auxiliary  saquenca 

P 

c(n)  •  r  '*'(n)  I  a.  r  (n-k)  (23) 

*  k-1  * 

in  which  the  aj^  coefficients  as  generated  from 
expression  (19)  are  used.  Due  to  the  nature  of  the 
causal  image  sequence  and  the  underlying  Yule- 
Walker  equations  (17) ,  it  is  readily  shown  that 
this  auxiliary  sequence  is  identically  zero  outside 
the  indexing  range  0  ^  n  ^  max(p,q).  Using  this 
fact,  the  Fourier  transform  of  relationship  (23) 
is  found  CO  yield 
s 

C  (e^“)  •  I  c(n)e"^“"  ,  s  »  Bax(p,q)  (24) 

*  n"0 

•  A  (eJ“)S  ^(.J")  (25) 

p  X 

Using  this  result  in  equation  (22),  Che  required 
AXfH  spectral  density  formulation  Is  obtained 


,  A  (eJ")c\e>)  +A*(eJ")C.(.J") 

c  fa-*")  --2 _ i - e - S - 

*  A  (*J“)A*(.J“) 

P  P 


To  obtain  the  ARHA  model's  bj^  coefficients, 
we  next  Incorporate  relationships  (4)  and  (26)  to 
generate  the  complex  conjugate  sysaecrlcal  poly¬ 
nomial  (In  the  eJ“)  expression 


B  (eJ")B*(.>)  -  A  (e>)C^.J“)  A*(e>)C  (.>) 

a  g  0  8  p  8 


(27) 


A  spsccral  factorazatlon  of  the  right  side  poly¬ 
nomial  will  yield  2q  roots  which  occur  in  complex 
conjugate  reciprocal  sets.  One  then  need  only 
salacc  an  appropriate  q  of  Chase  roots  to  detar- 
mlna  the  required  Bq(eJu)  term  (e.g.,  tha  minimum 
phase  selection) ■ 


To  summarize,  the  spectral  density  function 
and  tha  modal  coefficients  corresponding  to  an 
ARHA  time  series  of  order  (p,q)  may  ba  obtained  by 
following  the  systematic  procedure  outlined  In 
Table  1.  To  carry  out  this  process.  It  is 
necessary  to  have  knowledge  of  the  order  pair 
(p,q)  and  the  autocorrelation  elamoncs  r^Cn)  for 
0  <  n  i  q+p. 
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1. 

Solve  Relationship  (19)  for  Che  p  auto¬ 
regressive  ai[  coefficients.  This  will  re¬ 
quire  setting  c  ^  p. 

2. 

Generate  Che  auxiliary  sequence  c(n)  and  its 
Fourier  transform  using  expressions  (23)  and 

(24),  respectively. 

3. 

The  required  spectral  density  is  given  by 
relationship  ( 26) . 

4. 

Perform  a  spectral  factorization  of  Che 
polynomial  B(eJo))B*(eJ'^)  as  given  by  ex¬ 
pression  (27)  to  obtain  the  required  b^ 
coefficients. 

Table  1:  Generation  of  the  spectral  density  and 
Che  ARKA  model  parameters  associated  with  a 
given  set  of  autocorrelation  values. 

V.  EXAMPLE 

To  desionstrate  the  relative  effectlvenesa  of 
Che  MA,  AR,  and  ARMA  modeling  schemes  presented  in 
Che  lasc  chree  sections,  let  us  consider  the  class¬ 
ical  problem  of  generating  spectral  estimates  of 
a  time  series  composed  of  two  sinusoids  in  additive 
white  noise.  Namely,  the  time  series  will  be 
governed  by  Che  relationship 

x(n)  «  aj^3in(<>ij^n+dj^)  +a2Sla(<u2tt+92)  +w(n)  (28) 

in  which  and  $2  arc  independent  random  variables 
uniformly  distributed  on  [0,2n]  and  v(n)  Is  a  zero 
mean  white  noise  process  whose  individual  elemeots 
have  variance  The  sinusoidal  amplitudes  a^  and 
a2,  and,  normalized  frequencies  wi  and  U2  will  be 
here  Caken  co  be  unknown  constants.  ^ 


It  is  readily  shown  chat  the  autocorrelation 
sequence  corresponding  to  this  time  series  is 
given  by 

,2  .2  2 

r  (n)  •  1  co8(u,n) -f  2  cos(u,n) +o^S(n)  (29) 

X  ^  1  -y-  2 

Upon  caking  Che  Fourier  transform  of  this  suto- 
correlation  sequence,  the  associated  spectral  den¬ 
sity  function  is  found  to  be 
(  a  2 

»  1  nf <(M-M^^^.d(urHi2l  1 

^a  2  2 

>  _2_n  (d ((tf— )+d (urh*)-)  J  +0 
2  2  1 

for  -il  (30) 

This  spectral  density  function  is  seen  to  be  com¬ 
posed  of  dirac  delta  functions  located  at  fre¬ 
quencies  t  Lj;  &  f  UI2  riding  on  Cop  of  a  constant 
value  due  Co  Che  additive  white  noise. 


choices  for  the  time  series  parazMters  were  caken 
to  be  aj^  “  >^20,  uj  •  0,01,  a2  “  fT,  <J2  "  0.42611  and 
o2  a  1.  This  selection  provides  individual  sinu¬ 
soid  slgnsl-co-noise  ratios  of  lOdB  (decibels)  and 
OdB.  Due  to  Che  relative  closeness  of  the  sinu¬ 
soid  frequencies,  this  example  provides  an  excel¬ 
lent  measure  of  the  frequency  resolution  capabili¬ 
ties  of  Che  three  modeling  procedures.  A  brief 
description  of  Che  results  obtained  for  this 
exaaqile  now  follows. 

MA  Spectral  Estimstea 

The  autocorrelation  elements  as  specified  by 
expression  (29)  where  incorporated  into  the  HA 
spectral  model  relationship  (6)  in  which  Che 
window  function  is  caken  co  be  rectangular  (l.e. , 
w(n)  •  1  for  0  £  n  £  q) .  The  spectral  estimates 
thereby  achieved  for  the  specific  order  selections 
q  a  IS,  30,  and  200  are  displayed  in  Figure  1. 

^om  these  plots  it  la  clear  chat  the  MA  modeling 
procedure  la  unable  co  resolve  the  sinusoids  for 
Che  order  selections  q  •  is  and  30.  Whan  q  is  set 
CO  200,  it  is  possible  co  just  barely  detect  Che 
presence  of  the  lower  amplitude  sinusoid.  It  is 
apparent  from  this  example  chat  classical  Fourier 
approaches  provide  relatively  poor  vehicles  for 
achieving  frequency  resolution  even  when  exact 
autocorrelation  elemenca  are  used. 

AR  Spectral  Estimates 

The  autocorrelation  elements  (29)  were  next 
incorporated  into  the  optimum  one-step  predictor 
(or  itaxlssim  entropy)  expression  (IS) ,  with  AR 
order  choices  of  p  •  IS  and  30.  nie  evo  AR  spect¬ 
ral  esclmatas  which  resulted  are  shown  in  Figure 
2  where  it  is  apparent  chat  a  frequency  resolution 
is  achieved  for  p  <•  30,  but,  not  for  p  •  IS.  In 
contrast  to  the  classical  Fourier  approach,  Che 
AR  spectral  modeling  procedure  is  capeble  of 
achieving  the  required  frequency  resolution  with 
a  reasonably  small  order  model. 

ARMA  Spectral  Estimate 

In  Che  final  modeling  approach,  the  auto¬ 
correlation  elements  (29)  were  next  used  co  obtain 
an  ARMA  modal  of  order  (4,4)  using  the  procedure 
as  outlined  in  Table  1.  The  resultant  spectral 
escimaca  is  plotted  in  Figure  3  and  is  seen  co 
correspond  precisely  with  chat  as  given  in 
equation  (30). 3  This  should  not  be  surprising 
considering  the  fact  that  the  given  autocorrelation 
sequence  (29)  is  an  ARMA  time  series  of  order  (4,4). 
Thus,  the  model  used  in  this  case  precisely  mstchas 
the  time  series  being  examined. 

VI.  HIGH  PERFORMANCE  METHOD  OF 
ARMA  SPECTRAL  MODELING 


Using  Che  given  autocorrelation  elamanCs  (29) 
as  entries,  the  Chree  spectral  asclmatlon  pro¬ 
cedures  Just  described  were  next  utilized  to 
generate  MA.  AR,  and  ARMA  models.  The  specific 


^The  normalized  frequency  variables  are  caken  co 
lie  In  Che  inteival  [0,n]. 


It  is  possible  to  adapt  many  of  the  ideas  of 
Section  IV  to  achieve  an  ARMA  spectral  estimate 
whan  only  the  time  series  observation  (3)  are 


If  Infinite  precision  computations  ware  utilised, 
there  would  be  an  exact  pole-zero  cancellation  and 
an  associated  constant  spectral  density  function. 
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available  (and  not  autocorrelation  values) .  He 
shall  again  treat  separately  Che  cases  of  auto¬ 
correlation  and  moving  average  coefficient  deter¬ 
mination. 

Autoregressive  Coefficient  Eat<"*tlon 


ARMA  model  order  choice.  In  any  case,  the  system 
of  equations  with  these  estimate  subsclcuclona 
will  give  rise  to  Che  Cxl  Yule-Ualker  approximation 
error  vector  as  specified  by 

e  -  Y+X  a  +  Y+x  (36) 


To  Implemant  the  autoregressive  coefficient 
selection  process  as  represented  by  relationship 
(19) ,  It  will  be  necessary  to  cosqiute  appropriate 
autocorrelation  estimates  from  the  given  sat  of 
time  series'  observations.  The  high  performanca 
ARMA  method  effects  these  estimates  In  the  gulaa 
of  a  convenient  matrix  format  which  lends  Itself 
to  a  particularly  efficient  computational  realiza¬ 
tion  [2|&[4].  In  particular,  the  autocorrelation 
matrix  and  vector  required  In  expression  (19)  are 
estimated  according  to 

R  •  Y^X 

(31) 

i.  * 

(32) 

where  the  dagger  symbol  denotes  the  operation  of 
complex  conjugate  transposition.  The  (N-p)«p 
Toeplitz  type  matrix  X  la  specified  by 

'xlp)  X(p-l)4,  .  . 

x(l)  1 

x(p-t-l)  x(p)  .  .  . 

X  -  1  . 

x(2)  1 

1 

.  1 

1 

(33) 

'X(N-L)  x(N-2) 

.  i 

x(N-p)  i 

J 

while  the  (N-p)«t  Toeplitz 
form 

fx(p-q)  x(p-q-l)  . 

type  matrix  Y 

.  .x(p-q-C+-l)‘| 

has  the 

.x(p-q-*-l/  x(p-q)  .  . 

.  .x(p-q-t+2)! 

1 

Y  • ! 

1  (34) 

jx(N-q-l)  x(N-q-2)  . 

.  .x(N-q-c)  ! 

and  X  is  a  (N-p)-"!  vector 

given  by^ 

X  -  lx(p+-l) ,  x(p+2) , 

.  .  .  x(N)l’ 

(35) 

In  formulating  matrix  Y,  we  have  used  the  con¬ 
vention  of  setting  to  zero  any  elements  x(k)  for 
which  k  lies  outside  Che  observation  Index  range 
1  ^  k  <  N. 

If  Che  autocorrelation  matrix  and  vector 
esclmaces  (31)  and  (32),  respectively,  are  substi¬ 
tuted  Into  the  Yule-Walker  relationship  (19), 
however.  It  Is  generally  found  chat  the  resultant 
system  of  c  equations  In  the  p  autoregressive 
coefficients  Is  Inconsistent  for  c  '  p.  This  la 
due  CO  Inevitable  Inaccuracies  In  the  auto¬ 
correlation  estimates,  and,  to  a  possible  Improper 
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A  more  generalized  version  of  this  aaclmstlon 
scheme  can  be  obtained  by  substituting  the  Inte¬ 
ger  k  for  p  wherever  p  appears  In  relationships 
(33)  -  (35).  For  ease  of  presentation,  k  la 

here  restricted  to  be  p. 


Upon  caking  the  expected  value  of  £,  It  Is  found 
chat  for  the  ASHA  modeling  order  choice  In  which 
q  ^  p,  chat  this  expectation  results  In 


E{e(k)} 


(N-q 


p  P 

-k)  r^(q+k)  +  [ 


m>l 


l£k£t 

(37) 


idtUe  for  Che  modeling  order  case  q  <  p  this  ex¬ 
pectation  produces 


E(e(k)>  - 


r 

r  (q-fk)  +  y  a_r 


(N-p)  r  (q-fk)  +  [  a  r  (q^-k-m) 


m*l 


(N-q-k) ir  (q+k)  +  ^a  r  (q+k-m) 


,  l£kj^p-q 

,  p-q  <  k  <_t 
(38) 


In  either  ordering  case,  It  Is  seen  that  when  the 
time  series  la  an  ASMA  process  of  order  (p,q),  the 
expected  value  of  the  error  vector  £  can  be  made 
equal  CO  zero  by  a  proper  choice  of  the  auto¬ 
regressive  coefficient  vector  a  .  Namely,  this 
selection  would  be  such  that  the  underlying  Yule- 
Walker  equations  (17)  are  satisfied.^  This  Implies 
chat  the  system  of  equations  (36)  with  £,  •  pro¬ 
vides  an  unbiased  and  a  consistent  ssclmaCe~of  the 
Yula-Walkcr  equations. 


With  the  above  thoughts  In  lalnd,  an  appealing 
approach  to  selecting  the  autoregressive  co¬ 
efficient  vector  is  iMMdlataly  suggested.  Namely, 
a  is  chosen  so  as  to  make  the  error  vector  "as 
close"  to  Its  expected  value  of  £  as  possible. 

This  Is  of  course  predicated  on  the  assumption 
that  the  time  series  Is  an  ARMA  process  of  order 
(p,q)  or  less.  In  order  to  attain  a  tractable 
procedure  for  selecting  an  appropriate  auto¬ 
regressive  coefficient  vector,  we  shall  Introduce 
the  following  quadratic  functional 

f(a)  -  e^A  s  (39) 


In  which  A  Is  a  txt  positive-definite  diagonal 
matrix  with  diagonal  elements  that  Is 
Introduced  In  order  to  provide  one  with  the  option 
of  weighting  differently  the  various  error  vector 
components.  It  Is  a  simple  matter  to  show  that  an 
autoregressive  coefficient  vector  which  will 
render  this  quadratic  functional  a  mlnlsnim  must 
satisfy 

X^Y  AY+X  a*  •  -  X*Y  AY+x  (40) 


1 -  i 

A  little  though  will  convince  oneself  that  this 

same  conclusion  will  be  reached  If  both  q  and  p  j 

are  at  least  equal  to  tha  numerator  and  danoml-  j 

nator  orders,  respectively,  of  tha  underlying  j 

ARMA  time  series.  I 
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On*  Chan  slaply  solve*  eUs  consiscant  systaa  of 
p  linaar  aquations  in  tha  p  unknovn  autoragcaaalve 
coafflclancs  co  obtain  an  astlaac*  for  tha  danoai- 
nacor  dynaaica  of  ch*  ARMh  loodal. 

Movina  -  Avaraaa  Coafflcianc  Eaciaacion 

Thara  exist  several  procedures  for  astlaacing 
ch*  ASMA  aodal's  moving  average  coafflcients.  Ua 
shall  now  briefly  describe  two  procedures  which 
have  provided  satisfactory  perforaance  and  in  a 
sense  coapleaant  one  another. 

(i)  Cj^  Method 

The  procedure  which  has  provided  the  beat 
frequency  resolution  behavior  is  a  direct 
adaption  of  the  C|^  aethod  as  described  in 
Section  IV  (and  ref.  [3])‘  In  particular, 
using  the  sec  of  aucoregreaslve  coefficient 
estlaacas  as  obtained  froa  expression  (40) 
and  a  suitable  set  of  autocorrelation 
esclaatea  r((o)  for  n  •  0,1, . . . ,Bax(q,p)  one 
coapuces  the  c(a)  coefficients  using  expres¬ 
sion  (23) .  These  coefficients  are  then  used 
to  achieve  tha  desired  ARMA  spectral  estlaace 
whan  incorporated  into  ralatlonahlp  (24)  and 
ultlaataly  relationship  (26).  Although  pro¬ 
viding  an  excellent  frequency  resolution 
behavior,  this  procedure  suffers  ch*  drawback 
of  not  having  a  guaranteed  oonnegaciv*  defin¬ 
ite  spectral  density  function.^  It  is  with 
this  In  Bind  chat  the  following  well  known 
saooched  perlodograa  aethod  was  adapted  [13]. 


Cj^(n)  •  £(a +p  1 -t-kd)  O^a^q  (43) 

0  <k<L-l 

where  "d"  is  a  positive  integer  which  specifies  the 
tlae  shift  between  adjacent  sagaants.  These  indi¬ 
vidual  segaents  will  overlap  if  d  ^q  and  will  per¬ 
fectly  partition  the  residual  sequence  when  d  ■  q  1. 
In  order  to  include  only  coaputed  eleaancs,  tha 
relevant  paraaeters  must  be  selected  so  that 
q  4.p  4-1  +  (L-l)d  Next  the  perlodograa  for  each 

of  these  L  sagaants  is  taken  and  then  averaged  to 
obtain  the  desired  qfh  order  saooched  perlodograa, 
chat  is 

S  (eJ")  I  w(n)c,.(n)e--'““I  f  (44) 

®  ^  k-0 


q 

2 

1 

q+1 

[  w(n)E.  (n)e“^“” 
n-0 

where  w(n)  Is  a  window  sequence  chat  is  normally 
selected  to  be  rectangular  (l.e.,  w(n)  •  l//q+l  for 
O^n^q).  It  is  readily  shown  chat  this  procedure 
results  in  a  desired  nonnegatlve  qfh  order  HA 
spectral  density  estimate.  Unfortunately,  its 
frequency  resolution  capability  is  generally  not 
of  tha  same  quality  as  chat  of  tha  c^  aethod. ^ 

On  ch*  other  hand,  Che  saooched  parlodograa  aethod 
provides  more  saoothly  behaved  spectral  estlaataa 
which  contain  fewer  spurious  effects. 

To  auaaarlz*,  tha  required  ASMA  spectral  aodal 
is  obtained  by  following  the  syataaatlc  procedure 
outlined  in  Table  2.  The  nuaerator  dynamic  estl- 
aatlon  procadura  co  be  used  will  of  course  depend 
on  the  particular  characteristic  being  sought  (e.g., 
frequency  resolution,  smoothness,  etc.) 


(11)  Smoothed  Perlodograa  Method 


In  Che  aaoothed  perlodograa  approach, 
one  first  coapuCas  the  so-called  "residual” 
clae-serles  aleaent*  according  to  the  relatlon- 
shlp  (see  ref.  (3]). 

P  ^  . 

£(n)  ~x(n)  +  ^  a,  x(a-k)  for  p<  n<  H  (41) 

k-1  “ 


in  which  ch*  iyT  aucoregreaslve  coefficients 
as  obtained  by  solving  expression  (40)  are 
Incorporated.  Froa  this  relationship  the 
spectral  density  expression  directly  follows 


S,(.>) 


s,(.J“) 


(42) 


If  S,(*J“)  is  to  correspond  co  an  A8HA  spect¬ 
ral  modal  of  order  (p,q),  it  is  clear  that  a 
qth  order  NA  spectral  astlaata  for  the  resi¬ 
dual  spectral  density  S((*J<^)  auat  be  obtained 
and  Chen  substituted  into  relationship  (42). 
Tha  smoothed  perlodograa  has  been  found  co  be 
a  useful  Cool  for  this  purpose. 


In  the  smooched  parlodograa  method,  one  first 
partitions  the  computed  residual  alsaanta  (41)  into 
L  secants  each  of  length  q^l  as  specified  by 


1. 

Specify  values  for  the  ARMA  model's  order 
parameter  pair  (q,p),  the  Tule-Walkcr 
equation  paraaater  t,  and,  Che  weighting 
aatrice' s diagonal  elements  Xuc* 

2. 

Using  the  tlae  series  observations  x(l), 
x(2) , . . . ,x(N) ,  construct  Che  aatrices  Z,  Y, 
and  vector  x  according  to  relationships  (33), 
(34),  and  (75),  respectively. 

3. 

Deceralning  the  model's  autoregressive  co¬ 
efficients  by  solving  relationship  (40) . 

4. 

The  nuaerator' s  dynaalcs  are  obtained  by 
using  either  the  (1)  cj^  aethod,  or, 

(11)  the  smoothed  parlodograa  aethod. 

Table  2:  Basic  steps  of  ch*  standard  high  per¬ 
formance  AINA  spectral  eatlaatlon  aethod: 

The  Block  Processing  Mod*. 

Tha  laproved  spectral  eaciaacion  perforaance 
obtained  in  using  this  high  performance  aethod 
over  contemporary  ARMA  techniques  such  as  the  Box- 
Janklns  aethod  is,  to  a  large  extent,  a  conse¬ 
quence  of  selecting  tha  integer  c  co  be  larger 
than  tha  alnlmsl  value  p.  With  Che  corresponding 
larger  set  of  Tula-Walker  equations  chat  are 


This  shorccoalng  asy  be  superficially  avoided  by 
taking  Che  absolute  value  of  ch*  spectral  estlaat*. 


^A  slallar  approach  share*  ch*  same  accrlbucas  as 
does  Ch*  smooched  perlodograa  [14], 
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chereby  being  approxlaaCed,  Ic  Intuitively  followa 
that  the  aodel's  autoregressive  coefficients  will 
be  less  sensitive  to  autocorrelation  estlaate 
errors  which  are  embodied  In  Y+X  and  Y'^'x  than  would 
be  the  case  If  t  ware  sat  to  p  (as  In  the  Box- 
Jenklns  method) .  This  anticipated  Improvament  In 
spectral  estimation  behavior  when  using  the  high 
performance  method  has  In  fact  been  raallzed  on  a 
rather  large  number  of  numerical  examples  [ll-(4). 

It  is  shown  In  reference  [4]  i  (B)  that  this  high 
performance  method  also  lends  Itself  to  a  particu¬ 
lar  fast  adaptive  Implementation  mode  when  t  >  p . 
With  the  two  attributes  of  improved  spectral 
estimation  performance  and  computational  efficiency, 
this  new  procedure  promises  to  be  an  important 
spectral  estimation  tool. 

It  Is  of  Interest  to  note  that  when  q  -  0  and 
t  •  p,  Che  high  performance  ARtUi  spectral  esti¬ 
mation  method  reduces  to  the  well  luiown  AR  co- 
variance  method.  Moreover,  upon  letting  c  exceed 
p,  Che  resultant  sec  of  expanded  AR  Yule-Ualker 
equation  approximations  will  typically  result  in 
better  spectral  estimates  chan  the  standard  AR  co- 
variance  method.  To  the  auChorh  knowledge,  this  ap¬ 
proach  has  not  been  used  in  the  various  .AR  spectral 
estimation  procedures  developed  to  dace. 

VII.  CONCLUSION 

A  computationally  efficient  closed  form 
method  of  ARMA  spectral  estimation  has  been  pre¬ 
sented.  It  is  predicated  on  the  approximation  of 
a  set  of  Yule-Walker  equation  estimates  which  are 
generated  from  a  given  sat  of  time  series  obser¬ 
vations*  The  ARMA  model's  autoragrcasive 
coefficients  are  determined  by  solving  a  consistent 
system  of  linear  equations. 

The  spectral  estimation  performance  of  this 
ARMA  modeling  procedure  has  been  empirically  found 
CO  exceed  chat  of  such  counterparts  as  Che  maximum 
entropy  and  Box-Jcnklns  methods  (e.g.,  see  refs. 
[1]-14]  &  19)).  This  behavior  Is,  to  a  large 
extent,  a  consequence  of  the  fact  that  more  chan 
Che  minimal  number  of  Yule-Walker  equation  esti¬ 
mates  are  being  approxlsiated  to  obtain  the  result¬ 
ant  ARMA  model  parameters. 
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